#include "mhd.hh"
#include <math.h>
#include <dune/constant-unit/constant.hh>
namespace Dune {
  namespace Initial {
    MHD::MHD(Dune::ParameterTree& ptree){
      set_initcond_para("dim",1.0);
      set_initcond_para("nvar",1.0);
      set_initcond_para("sod_shocktube_leve",0.0);
      set_initcond_para("sod_shocktube_stone",0.0);
      set_initcond_para("twoshock_leve",0.0);
      set_initcond_para("briowu",0.0);
      set_initcond_para("rj2a",0.0);
      set_initcond_para("rj4d",0.0);
      set_initcond_para("blastwave",0.0);
      //
      set_initcond_para("cen1",0.0);
      set_initcond_para("rhoI",0.0);
      set_initcond_para("rhoII",0.0);
      set_initcond_para("mom1I",0.0);
      set_initcond_para("mom1II",0.0);
      set_initcond_para("mom2I",0.0);
      set_initcond_para("mom2II",0.0);
      set_initcond_para("mom3I",0.0);
      set_initcond_para("mom3II",0.0);
      set_initcond_para("b2I",0.0);
      set_initcond_para("b2II",0.0);
      set_initcond_para("b3I",0.0);
      set_initcond_para("b3II",0.0);      
      set_initcond_para("eI",0.0);
      set_initcond_para("eII",0.0);
      set_initcond_para("gamma",5.0/3.0);
      set_initcond_para("b1zero",0.0);  
      set_initcond_paras(ptree);
    }

    MHD::~MHD(){
    }

    void MHD::set_initcond_paras(Dune::ParameterTree& ptree){
      para["sod_shocktube_leve"]=ptree.get("initial.sod_shocktube_leve",0.0);
      para["sod_shocktube_stone"]=ptree.get("initial.sod_shocktube_stone",0.0);
      para["twoshock_leve"]=ptree.get("initial.twoshock_leve",0.0);
      para["briowu"]=ptree.get("initial.briowu",0.0);
      para["rj2a"]=ptree.get("initial.rj2a",0.0);
      para["rj4d"]=ptree.get("initial.rj4d",0.0);
      para["blastwave"]=ptree.get("initial.blastwave",0.0);
      //
      para["dim"]=ptree.get("grid.dimension",0.0);
      para["cen1"]=ptree.get("initial.cen1",0.0);
      para["rhoI"]=ptree.get("initial.rhoI",0.0);
      para["rhoII"]=ptree.get("initial.rhoII",0.0);
      para["mom1I"]=ptree.get("initial.mom1I",0.0);
      para["mom1II"]=ptree.get("initial.mom1II",0.0);
      para["mom2I"]=ptree.get("initial.mom2I",0.0);
      para["mom2II"]=ptree.get("initial.mom2II",0.0);
      para["mom3I"]=ptree.get("initial.mom3I",0.0);
      para["mom3II"]=ptree.get("initial.mom3II",0.0);
      //
      para["b2I"]=ptree.get("initial.b2I",0.0);
      para["b2II"]=ptree.get("initial.b2II",0.0);
      para["b3I"]=ptree.get("initial.b3I",0.0);
      para["b3II"]=ptree.get("initial.b3II",0.0);
      para["eI"]=ptree.get("initial.eI",0.0);
      para["eII"]=ptree.get("initial.eII",0.0);
      para["gamma"]=ptree.get("initial.gamma",0.0);
      para["b1zero"]=ptree.get("initial.b1zero",0.0);
    }

    void MHD::initial_condition(double* location, double* holder){
      if(para["sod_shocktube_leve"]==1.0)
	sod_shocktube_leve(location, holder);
      if(para["sod_shocktube_stone"]==1.0)
	sod_shocktube_stone(location, holder);
      if(para["twoshock_leve"]==1.0)
	twoshock_leve(location, holder);
      if(para["briowu"]==1.0)
        briowu(location, holder);
      if(para["rj2a"]==1.0)
	rj2a(location, holder);
      if(para["rj4d"]==1.0)
	rj4d(location, holder);
      if(para["blastwave"]==1.0)
	blastwave(location, holder);
    }

    // Figure 14.3 of book: Leveque R. J.
    // similar with Fig 8 of Stones J. M. 2008ApJ
    void MHD::sod_shocktube_leve(double* location, double* holder){
      para["cen1"] =0.5;
      holder[0] = (location[0]<para["cen1"]) ? 3.0 : 1.0;
      holder[1] = 0.0;
      holder[2] = 0.0;
      holder[3] = 0.0;
      holder[4] = (location[0]<para["cen1"]) ? 
        3.0/(para["gamma"]-1.0) : 1.0/(para["gamma"]-1.0);
      holder[5] = 0.0;
      holder[6] = 0.0;
      return;
    } 

    // Fig 8 of Stones J. M. 2008ApJ
    void MHD::sod_shocktube_stone(double* location, double* holder){  
      para["cen1"] =0.5;  
      para["gamma"]=1.4;
      holder[0] = (location[0]<para["cen1"]) ? 1.0 : 0.125;
      holder[1] = 0.0;
      holder[2] = 0.0;
      holder[3] = 0.0;
      holder[4] = (location[0]<para["cen1"]) ? 1.0/(para["gamma"]-1.0) : 0.1/(para["gamma"]-1.0);
      holder[5] = 0.0;
      holder[6] = 0.0;
      return;
    }

    // Figure 14.4 of book: Leveque R. J.
    void MHD::twoshock_leve(double* location, double* holder){
      para["cen1"] =0.5;  
      para["b1zero"]=0.0;
      holder[0] = (location[0]<para["cen1"]) ? 1.0 : 2.0;
      holder[1] = (location[0]<para["cen1"]) ? 1.0*3.0 : 2.0*1.0;
      holder[2] = 0.0;
      holder[3] = 0.0;
      holder[4] = (location[0]<para["cen1"]) ? 4.5+1.0/(para["gamma"]-1.0) : 1.0+1.0/(para["gamma"]-1.0);
      holder[5] = 0.0;
      holder[6] = 0.0;
      return;   
    }

    // Figure 13 of Stones J. M. 2008ApJ
    void MHD::briowu(double* location, double* holder){
      para["cen1"] =0.5;  
      para["b1zero"]=0.75;
      para["gamma"]=2.0;
      holder[0] = (location[0]<para["cen1"]) ? 1.0 : 0.125;
      holder[1] = 0.0;
      holder[2] = 0.0;
      holder[3] = 0.0;
      holder[4] = (location[0]<para["cen1"]) ? 
        0.5*(1.0+std::pow(para["b1zero"],2))+1.0/(para["gamma"]-1.0) : 
        0.5*(1.0+std::pow(para["b1zero"],2))+0.1/(para["gamma"]-1.0);
      holder[5] = (location[0]<para["cen1"]) ? 1.0 : -1.0;
      holder[6] = 0.0;
      return;
    }

    // Figure 14 of Stones J. M. 2008ApJ
    // RJ2a means fig2a of Ryu Dongsu & Jones TW 1995ApJ
    void MHD::rj2a(double* location, double* holder){
      para["cen1"] =0.5;    
      para["b1zero"]=2.0/std::sqrt(4.0*Dune::Constant::dpi);
      para["gamma"]=5.0/3.0;
      holder[0] = (location[0]<para["cen1"]) ? 1.08 : 1.0;
      holder[1] = (location[0]<para["cen1"]) ? 1.08*1.2 : 0.0;
      holder[2] = (location[0]<para["cen1"]) ? 1.08*0.01 : 0.0;;
      holder[3] = (location[0]<para["cen1"]) ? 1.08*0.5 : 0.0;;
      holder[5] = (location[0]<para["cen1"]) ? 
        3.6/std::sqrt(4.0*Dune::Constant::dpi) : 4.0/std::sqrt(4.0*Dune::Constant::dpi);
      holder[6] = (location[0]<para["cen1"]) ? 
        2.0/std::sqrt(4.0*Dune::Constant::dpi) : 2.0/std::sqrt(4.0*Dune::Constant::dpi);
      holder[4] = (location[0]<para["cen1"]) ?
        0.5*(std::pow(para["b1zero"],2)+3.6*3.6/(4.0*Dune::Constant::dpi)+
        2.0*2.0/(4.0*Dune::Constant::dpi))+0.5*1.08*(1.2*1.2+0.01*0.01+0.5*0.5)+
        0.95/(para["gamma"]-1.0) :
        0.5*(std::pow(para["b1zero"],2)+4.0*4.0/(4.0*Dune::Constant::dpi)+
        2.0*2.0/(4.0*Dune::Constant::dpi))+1.0/(para["gamma"]-1.0);
      return; 
    }
    
    // RJ4d: Figure 15 of Stones J. M. 2008ApJ
    void MHD::rj4d(double* location, double* holder){
      para["cen1"] =0.5;  
      para["b1zero"]=0.7;
      para["gamma"]=5.0/3.0;
      holder[0] = (location[0]<para["cen1"]) ? 1.0 : 0.3;
      holder[1] = (location[0]<para["cen1"]) ? 0.0 : 0.0;
      holder[2] = (location[0]<para["cen1"]) ? 0.0 : 0.0;
      holder[3] = (location[0]<para["cen1"]) ? 0.0 : 0.3*1.0;
      holder[5] = (location[0]<para["cen1"]) ? 0.0 : 1.0;
      holder[6] = (location[0]<para["cen1"]) ? 0.0 : 0.0;
      holder[4] = (location[0]<para["cen1"]) ?
        0.5*(std::pow(para["b1zero"],2))+1.0/(para["gamma"]-1.0) :
        0.5*(std::pow(para["b1zero"],2)+1.0)+0.5*0.3+0.2/(para["gamma"]-1.0);
      return;
    }

    // blastwave for 2D/3D MHD test
    void MHD::blastwave(double* location, double* holder){
      int dim = (int)para["dim"];
      para["rhoI"]=1.0;
      para["pI"]=3.0;
      para["pII"]=1.0;
      para["radius"]=0.15;
      double bzero=1.0/1.415926;
      double distance=0.0;
      double benergy=0.5;
      holder[0] = para["rhoI"];
      holder[1] =0.0;
      holder[2] =0.0;
      holder[3] =0.0;
      for(int i=0; i<dim; i++)
        distance+=std::pow(location[i]-0.5,2);
      distance=std::sqrt(distance);
      holder[4] = distance<=para["radius"] ? para["pI"]/(para["gamma"]-1.0)+benergy :
        para["pII"]/(para["gamma"]-1.0)+benergy;
      holder[5] =bzero;
      holder[6] =bzero;
      holder[7] =0.0;
      return;
    }
 } // Initial
}  // Dune
